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ABSTRAGT 

In this work, we present, for the first time, a numerical study of the Bondi-Hoyle accretion with 
density gradients in the fully relativistic regime. In this context, we consider accretion onto a Kerr 
Black Hole (BH) of a supersonic ideal gas, which has density gradients perpendicular to the relative 
motion. The set of parameters of interest in this study are the Mach number, A^, the spin of the BH, 
a, and the density-gradient parameter of the gas, e^. We show that, unlike in the Newtonian case, 
all the studied cases, especially those with density gradient, approach a stationary flow pattern. To 
illustrate that the system reaches steady state we calculate the mass and angular momentum accretion 
rates on a spherical surface located almost at the event horizon. In the particular case of Ai = 1, 
Cp = 0.5 and BH spin a = 0.5, we observe a disk-like configuration surrounding the BH. Finally, we 
present the gas morphology and some of its properties. 

Subject headings: black hole physics - accretion - relativistic processes - hydrodynamics - meth¬ 
ods: numerical 


1. INTRODUGTION 

1.1. Accretion Onto Astrophysical Black Holes 

Stellar-Mass Black Holes (sBHs): Astrophysical black 
holes (BHs) come in a variety of masses and spins. sBH 
candidates with AMq < Mbh ^ SOMq have been ob¬ 
served (see, e.g., table 1 in Moreno Mendez 2013, and ref¬ 
erences therein). The masses of the BHs can be obtained 
through two main channels. The first of which is by col¬ 
lapse of the progenitor star and capture of fallback mass. 
The second channel is through mass transfer from a com¬ 
panion star. Of course, other channels are possible; e.g., 
the BH could be the merger of two compact objects (NS- 
NS, BH-NS, WD-WD, NS-WD, BH-WD, etc.), a com¬ 
pact object and a main sequence (MS) star or through a 
common-envelope phase; the last two mechanisms could 
produce a Thorne-Zytkow Objects (TZO) that eventu¬ 
ally forms a BH (Thorne & Zytkow 1975). 

Besides the masses of sBHs it has also been possible 
to obtain their spins. These may reveal important fea¬ 
tures of the BHs history. Some of these spins^ have been 
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^ The spin, or Kerr parameter, is defined as a* = m%h ~ 

(where is the angular momentum of the BH and Mbh 

is its mass) and it varies between —1, for accretion disks which 
counterrotate with respect to the BH, to 0 for nonspinning BHs, 
to +1 for corrotating accretion disks. 


estimated through three different methods and seem to 
be distributed between ^ —0.2 and > 0.98. All these 
methods depend on observing phenomena of an accretion 
disk surrounding the sBH, thus a donor star is neces¬ 
sary. One method for measuring sBH spins is by fitting 
the X-ray continuum, (see McGlintock et al. 2015, for 
a review). A second method is through X-ray reflection 
spectroscopy (of the Fe K-a line; see Reynolds 2015, for a 
review). The third method that has been used is through 
Quasi-Periodic Oscillations (QPOs) from the precession 
of the accretion disk (see, e.g., Axelsson et al. 2005). 

Stellar-mass BHs may acquire their spins through a 
variety of mechanisms (Lee et al. 2002; Yoon & Langer 
2005; Woosley & Heger 2006). In most BHs known in 
Low-Mass X-ray Binaries (LMXBs), the spins may be the 
result of pre-BH formation processes, either by accretion 
through mass transfer (Moreno Mendez et al. 2011), or 
through tidal-synchronization after a post-Gase-C-mass- 
transfer-and-common-envelope episode (Brown et al. 
2007; Moreno Mendez 2014); this mechanism may even 
lead to the production of gamma-ray burst and hyper¬ 
nova events (e.g.. Brown et al. 2008; Moreno Mendez 
et al. 2014; but see also Fragos & McGlintock 2014). 

For the BHs in High-Mass X-ray Binaries (HMXBs), 
the spin estimates seem to pile on the high end, i.e., 
they vary between 0.84 and > 0.98. Stellar evolution 
in binaries has trouble explaining said spins and so do 
core-collapse mechanisms (Moreno Mendez & Gantiello, 
in progress), thus Moreno Mendez et al. (2008); Moreno 
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Mendez (2011) have suggested that wind-driven mass 
transfer with hypercritical accretion may be necessary 
to explain these binaries. 

Intermediate-Mass Blaek Holes (IMBHs): IMBHs are 
the logical intermediary between sBHs and supermassive 
BHs (SMBHs). If the later are produced by accretion 
onto the former or from mergers of massive stars or sBHs, 
then IMBHs should be abundant. On the other hand, if 
they are only produced from core-collapse of Population 
III stars, they may have already turned into SMBHs. 
Indeed, they are excellent candidates to explain the ob¬ 
served Ultra Luminous X-ray (ULX) sources. Nonethe¬ 
less, many ULXs have been identified with LMXBs (Ba- 
chetti et al. 2014) and HMXBs (Motch et al. 2014; Liu 
et al. 2013). There seems to be, however, a good candi¬ 
date of an IMBH in a ULX where, using High-frequency 
QPOs (with 3:2 ratio), Pasham et al. (2014) use inverse¬ 
scaling of sBH as well as a relativistic precession model 
to determine a mass of Mbh — 4OOM0. 

Super-Massive Blaek Holes (SMBHs): Explaining the 
mass of supermassive black holes (SMBHs) is one of the 
most interesting current problems of astrophysics. A 
standard approach to the problem assumes that these 
holes are the result of the accretion onto intermediate- 
mass black-hole (IMBH) seeds. This raises the question 
of where such seeds come from in the first place. The 
answer to this question has opened a wide set of possibil¬ 
ities. It could be that disks formed in protogalaxies allow 
the infall of matter that collapses to form black holes of 
10 ^Mq (Koushiappas et al. 2004). Alternatively, seeds 
of 10^ Mq can be formed due to the runaway collision of 
compact stellar clusters in low metallicity protogalaxies 
at z ^ 10 — 20 that additionally allow accretion during 
the quasar era (Devecchi & Volonteri 2009). Also, seeds 
can be the result of the core collapse of dense clusters and 
form 10 ^Mq black holes (Davis & Laor 2011). Another 
alternative is that seeds of 10^ Mq can be formed due to 
the collapse of very massive stars (Umeda et al. 2009) at 
extremely low metallicity. A more detailed model in this 
direction proposes that a supermassive primordial star 
forms in a region of the Universe with a high molecule- 
dissociating background-radiation field, and collapses di¬ 
rectly into a 10^ — IO^Mq black-hole seed (Johnson et al. 
2013). More recently, it has been proposed that primor¬ 
dial black holes, during the radiation-dominated era, can 
grow up to 10^ — IO^Mq (Lora-Clavijo et al. 2013b). 

Based on the study of the evolution of phase space dis¬ 
tributions, the growth process of seeds in standard anal¬ 
ysis of SMBHs considers that these are primarily fed by 
collisionless dark matter or stars (Lightman & Shapiro 
1977; Zhao et al. 2002). Previous results show that the 
timescale is extremely long for collisionless dark mat¬ 
ter to contribute significantly to black-hole growth (for 
instance Read & Gilmore 2003). In Guzman & Lora- 
Glavijo (201 lb, a) the conditions for stable accretion and 
runaway accretion of an ideal gas have been studied. 
Other more generalized cases, consider non-radial accre¬ 
tion processes which require the particles to overcome 
the angular momentum barrier in order to get the gas to 
the center of a galaxy (King & Pringle 2006). 

The assumption of self-interacting dark matter (SIDM) 
has been used to analyze the problem of SMBH growth. 
For instance in Balberg et al. (2002) and Balberg & 


Shapiro (2002) SIDM is used to model direct dynami¬ 
cal collapse and SMBH formation due to the gravother- 
mal catastrophe. Furthermore, this would explain the 
different SMBH seed masses in terms of the redshift at 
which the collapse took place. Munyaneza & Biermann 
(2005) show that fermionic dark matter can feed SMBH 
seeds to make them grow up to 10^ — 10 ^Mq. Instead, 
Saxton & Wu (2008) propose that the self-interaction 
is introduced using a polytropic equation of state for the 
pressure, thus, finding another mechanism for SMBH for¬ 
mation. Ostriker (2000) (and later Hu et al. 2006) ana¬ 
lyzed black hole formation due to the collapse of SIDM 
and also studied SMBH formation and growth due to 
the collapse and accretion of SIDM. Lora-Glavijo et al. 
(2014b) studied the accretion of a SIDM into a SMBH 
for the particular case of radial flows. They considered 
the evolution of space-time in order to have the formally 
correct black-hole growth rate. 

The characterization of SMBH spin is of vital impor¬ 
tance since it probes their growth history as well as 
their formation. In essence, scenarios in which SMBH 
growth is dominated by BH-BH mergers predict a popu¬ 
lation of modestly spinning SMBHs, whereas growth via 
gas accretion can lead to either, a rapidly-spinning, or 
a very slowly-spinning population (Moderski & Sikora 
1996; Volonteri et al. 2005; Reynolds 2015). SMBH spin 
can also be a potent energy source, and may well drive 
the powerful relativistic jets that are seen from many 
black hole systems (e.g., Blandford & Znajek 1977). 

1.2. Bondi-Hoyle-Lyttleton Aeeretion 

Bondi-Hoyle-Lyttleton (BHL) accretion deals with the 
evolution of a homogeneously distributed gas moving uni¬ 
formly toward a central compact object (Hoyle & Lyttle¬ 
ton 1939; Bondi & Hoyle 1944). Depending on whether 
the velocity of the gas is supersonic (or not) a shock 
cone is formed (or not). This process shows interesting 
properties when considered within the Newtonian and 
relativistic regimes, which have been explored based on 
several numerical studies. In the classical regime, which 
is ruled by Newtonian gravity, the most important sub¬ 
jects are the consequences on the morphology of the wind 
and the supersonic shocks that develop. A summary of 
results, under Newtonian gravity, can be found in (Edgar 
2004; Foglizzo et al. 2005). 

Unlike in the Newtonian regime, the relativistic ap¬ 
proach allows the study of BHL accretion in regions 
where the gravitational field is strong. Some studies 
in this direction have been carried out. The first one, 
performed by Petrich et al. (1989), studied the differ¬ 
ent accretion patterns developed by the relativistic gas 
during the accretion onto a BH. Later on, considering ax¬ 
ial and equatorial symmetries. Font & Ibanez (1998a, b); 
Font et al. (1998, 1999) reviewed the results obtained 
by Petrich et al. (1989) using more accurate methods. 
In the astrophysical context and using equatorial sym¬ 
metry, Donmez et al. (2011) showed that the shock cone 
vibrations can be associated with sources of QPOs. They 
also found a flip-flop type of unstable oscillation of the 
shock cone. However, it was later shown by Gruz-Osorio 
et al. (2012) that the flip-flop oscillation of the shock 
cone depends on the coordinates used to describe the ro¬ 
tating black hole, specifically, it was found that the flip- 
flop oscillation does not appear when Kerr-Schild coordi- 
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nates are used to describe the rotating black hole. More 
recently, considering axially symmetric fluxes, in Lora- 
Clavijo & Guzman (2013), the shock cone oscillations 
as a potential source of low and high frequency QPOs 
were studied. The ultrarelativistic BHL accretion onto a 
rotating black hole was recently reported considering ax- 
isymmetric fluxes (Penner 2013) and considering equato¬ 
rial symmetry (Cruz-Osorio et al. 2013). More realistic 
scenarios introduce astrophysically relevant ingredients 
like magnetic fields (Penner 2011), radiative terms (Zan- 
otti et al. 2011) and full general relativity in the context 
of supermassive black hole binary mergers (Farris et al. 
2010). 

This work uses penetrating coordinates (Kerr-Schild), 
which offer the possibily to place the excision inside the 
event horizon of the BH. It is well known that what hap¬ 
pens inside a BH, stays inside the BH, thus this further 
avoids the implementation of boundary conditions there. 

In the Newtonian regime, the BHL problem is attended 
by assuming the accretor is point-like and the size of the 
accretor is determined in terms of the accretion radius. 
In our case we are doing the exact opposite; this im¬ 
plies that if we want to resolve the accretor, i.e. the 
Schwarzschild radius must be similar in order of mag- 
nitud to the accretion radius, thus, the velocity of the 
wind must be high. This is so that the ratio of the ac¬ 
cretion radius (as well as that of the compact object) 
can be small enough as compared to the exterior numer¬ 
ical domain radius (the ratio must be at least an order 
of magnitude) so that no numerical artifacts affect the 
simulations Cruz-Osorio et al. (2012); Font et al. (1998). 
Hence, this allows to have good resolution and a reason¬ 
able computational timescale. The accretion radius is 
defined to help decide when a particle falls into the com¬ 
pact object. Such a scale is determined by the velocity of 
the wind and the equation of state of the gas. In our case 
we use the accretion radius only to choose the exterior 
numerical domain, that is about ten times the accretion 
radius. 

Another problem one may face is that the relativis¬ 
tic Euler equations may diverge or develop unphysical 
results; thus an atmosphere is implemented which pre¬ 
vents the rest mass density to be small enough such 
as to provoke said effects. The atmosphere we use is 
p = max{p^ 10“^^); such a value allows the convergence 
of our numerical methods. The need of the atmosphere 
is in fact one of the reasons why it is not -at the moment- 
trivial to simulate the evolution of fluids with ultralow 
density. Given these numerical limitations, the maxi¬ 
mum density difference between the BH^ and the non- 
homogeneous wind is of order Pbh/P w — 10^^. Also, the 
velocities in our models will be restricted to values be¬ 
tween 0.1c and 0.5c. This minimizes the possible scope 
of astrophysical scenarios that can be dealt with. 

1.3. Outlook 

In this manuscript we present for the first time a nu¬ 
merical model of relativistic BHL accretion onto a Kerr 

^ Here we understand by BH density the BH mass divided by the 
spherical volume using the radius of its event horizon, thus pbh oc 
paper we are working under the assumption that 
the space-time background is unaffected by the fluid contribution, 
i.e., we have a test fluid, which is consistent with our astrophysical 
models. 


BH of a non-homogeneous gas cloud. This problem was 
addressed initially, in the Newtonian regime, by Ruffert 
& Anzer (1995) and Ruffert (1997, 1999). They investi¬ 
gated the hydrodynamics of three-dimensional classical 
BHL accretion; especially the accretion of angular mo¬ 
mentum from a non-homogeneous medium and discussed 
some consequences for models of wind-fed X-ray sources. 
They found that the models with a density gradient ex¬ 
hibit non-stationary flow patterns, even though the Mach 
cone remains fairly stable. In more recent work, also in 
the Newtonian regime, MacLeod & Ramirez-Ruiz (2014) 
model the orbital inspiral of a neutron star (NS) through 
the envelope of its giant-branch companion during a com¬ 
mon envelope (GE) episode. They found that the pres¬ 
ence of a density gradient strongly limits the accretion 
by imposing a net angular momentum to the flow around 
the NS. 

Based on the limitations stated above, we accord¬ 
ingly list which models are appropriate for our differ¬ 
ent simulations which have been left unitless in order 
to maximize the posible scenarios in which they would 
be useful. Hence, in what follows, the mass of the 
BH shall be denoted M. Eollowing the usual conven¬ 
tion, we shall mostly use units in which G = c = 1. 
Thus, length and time are measured in units of M with 
1 M© = 1.48 X 10^ cm = 4.93 x 10“^ s. 

This paper is organized as follows: Section 2 proposes 
models for BHs of different masses. Section 3 describes 
the relativistic hydrodynamic equations as well as the nu¬ 
merical methods employed by our code. Section 4 shows 
our numerical results. And in section 5 we discuss them 
and produce our conclusions. 

2. MODELS 

We have produced a set of models which are BH mass 
independent, therefore, once the mass of the BH has been 
chosen, its density is established and one may obtain the 
wind density that hits the BH. The velocities of the flows 
hitting the BH are 0.1c, 0.2c, 0.3c, 0.4c and 0.5c. 

We next list a set of astrophysical situations where our 
model would become a good fit with the proper scaling. 

2.1. For Stellar-Mass BHs 

Relativistic winds and shockwaves: An astrophysical 
scenario where our models may apply for stellar-mass 
BHs is during a hypernova in a binary system. Suppose 
we have a short-orbital-period binary ( P < 0.4 day) 
with a massive WR star {Mwr > lOM©) and a ^ lOM© 
BH (e.g., Gyg X-3 may be a system in a situation similar 
to this description). The WR tidally synchronizes and 
rotates rapidly. Such a star will collapse into a compact 
object and will likely launch a GRB/HN explosion (see, 
e.g., Moreno Mendez et al. 2011). If ^ 5M0 collapse, 
the resulting BH could have > 0.8 and the available 
energy would be > 10^^ erg (or 1000 Bethes; 1 Bethe = 1 
B = 10^^ erg). If ^ 5M0 of material are expelled with 
a fraction (say 10% to 50% of this energy) their average 
velocity would be between 0.1c and 0.7c. At the time of 
the explosion, the orbital separation is about 3 P© to 4 
P©. Thus, the density of the ejecta at the orbit of the 
BH could be as large as one fifth of its density inside the 
star if the expansion is mostly equatorial or a hundredth 
if the expansion is spherically symmetric. This translates 
to densities as high as 10^ g cm“^ if material from the G 
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core is expelled. For a IOMq BH, the density is pbh ^ 
IQie g cm“^. Thus our model deals with py^ 10^ g 
cm“^, which is 1 to 2 orders of magnitude above what 
this model provides. Hence, our simulations give a good 
qualitative description of what should occur for such a 
scenario. 


2.2. IMBHs: 

IMBHs have a density range of pbh ^ 10^^ g cm“^ 
{Mbh = IO^Mq) to Pbh ^ 10^ g cm“^ {Mbh = 
10^Mq). Thus, given that our simulations are bound 
to p/ Pbh ^ 10“^^, we have a wind density regime that 
goes from 10 g cm“^ to 10“^ g cm“^. 

A scenario similar to the previous one would 
be ideal. I.e., an IMBH in the mass range of 
10^ Mq < Mbh ^ 10^ Mq, located, e.g., near the 
center of a globular cluster which were to tightly capture 
a massive WR star (maybe through a common envelope 
phase). If the orbital period was, again, 0.4 days, the 
orbital separation, d, would be 11 Rq ^ d < 50 R©, 
respectively. The WR radius being some 50 to 10 % of 
its Roche lobe (RL) radius (the BH radius being less 
than about a thousandth of its RL radius). For such 
an orbital separation the density of the ejecta should be 
around 5 x 10 “^ g cm“^ to 10 “^ g cm“^, at 11 Rq] 
and 10 “^ g cm“^ to 5 x 10 “^ g cm“^, at 50 Rq. Thus, 
these models are well within the range of our numerical 
calculations. 


2.3. SMBHs: 

SMBHs (IO^Mq to IO^Mq; although there is recent 
evidence for a SMBH as large as Mbh ^ 10^^ Mq, Lopez- 
Cruz et al. 2014) have a density range that goes from 
10^ g cm“^ to 10“^ g cm“^. The wind densities in our 
numerical simulations are, hence, confined between 10“^ 
g cm“^ and 10“^^ g cm“^. 

Among the possible scenarios where our models may 
be applied for SMBHs are, e.g., a SMBH binary where 
one of the BHs has relativistic jets pointing along the 
orbital plane. As the second BH steps into or out of the 
jet of the first BH, it is submerged in a stream of material 
with a density gradient. 

It is conceivable that the last model discussed for 
IMBHs (a nearby GRB pointing towards the BH) could 
similarly apply for SMBHs but at somewhat larger dis¬ 
tances. 

3. RELATIVISTIC HYDRODYNAMIC EQUATIONS AND 
NUMERICAL METHODS 

3.1. Relativistic Hydrodynamics Equations 

In order to solve numerically the relativistic Euler 
equations, we use the 3+1 decomposition of space-time, 
in which the space-time is foliated with a set of non¬ 
intersecting space-like hypersurfaces (see e.g. Alcu- 
bierre 2008; Baumgarte & Shapiro 2010; Rezzolla & Zan- 


otti 2013). The space-time is described with the line el¬ 
ement 


ds^ = —a^dt^ + ^ij{dx^ + P'^dt){dx^ + /^-^dt), (1) 

where a is the lapse function, are the shift vector 
components and ^ij are the components of the induced 
three metric that relates proper distances on the spatial 
hyper surfaces. 

The background space-time corresponds to a rotating 
black hole in Kerr-Schild coordinates, which allow us to 
place the inner boundary of the computational domain 
inside the event horizon. A discussion of the advantage 
on the use of these coordinates can be found in Cruz- 
Osorio et al. (2012). Once the geometrical elements of 
the space-time background are known, it is necessary to 
track the evolution of the fluid, for which we write down 
the general relativistic Euler equations. Eor a generic 
space-time these can be derived from the local conserva¬ 
tion of the stress-energy tensor 

V,(T^") = 0, (2) 

and the local conservation of the rest mass density 

v.(pw'^)=0, (3) 

where p is the proper rest mass density, is the four- 
velocity of the fluid and is the covariant derivative 
consistent with the four-metric of the space-time (1). 

We assume the matter field in the above equations is 
that of a perfect fluid with stress-energy tensor 

T^, = phu^u^^pg^\ (4) 

where p is the pressure, g^^ are the components of the 

four-metric and h the relativistic specific enthalpy given 
by h = 1 + e+ p/p, where e is the rest frame specific 
internal energy density of the fluid. 

It is well known that Euler’s equations develop discon¬ 
tinuities in the hydro dynamical variables even if smooth 
initial data are considered. Thus one may solve hydrody¬ 
namics equations using finite volume methods, as long as 
the system of equations is written in a flux balance law 
form, which in turn requires the definition of conservative 
variables. 

In order to obtain the general relativistic Euler equa¬ 
tions as a set of flux balance laws (Banyuls et al. 1997; 
Eont et al. 2000), we project the local conservation equa¬ 
tions along the space-like hypersurfaces and the normal 
direction to such hyper surfaces. A straightforward cal¬ 
culation yields the set of equations in the desired form 


1 


5t(\Aq) + di (v^f^*Hq)) = s(q), 


( 5 ) 


where g is the determinant of the four-metric (1), q is 
a vector of conservative variables, f^*^(q) are the fluxes 
along each spatial direction and s(q) is a source vector. 
These last quantities are given by: 
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q=[L), Mj, 
f«(q)= '' 


• /3* 

a 



[pr, phr'^Vj, phr'^-p-pTf, 


T 


T + v''p , 


(6) 

( 7 ) 


s(q) = [ 0 , - aT^'rO,]^. 


(8) 


In these expressions, 7 = det{^ij) is the determinant of 
the spatial metric, are the Christoffel symbols and 

v'^ is the 3-velocity measured by an Eulerian observer and 
defined in terms of the spatial part of the 4-velocity 
as v'^ = /T /^Y^, where E is the Lorentz factor given 

by r = 1/71 - 

It is still necessary to close the system of equations (5), 
for which an equation of state relating p = p(p, e) is used. 
We choose the gas to obey an ideal gas equation of state 

p = pe( 7 -l), (9) 

where 7 is the adiabatic index or the ratio of specific 
heats. Something to stand out is that the relativistic 
sound velocity Cg for an ideal equation of state can be 
written as = ^ 7(7 — 1 )/[P 7 — p (7 ~ 1 )], where its 
asymptotic value or its maximum permitted value is 
Cg = a /7 — 1- Thus, the choice of our initial values 
is restricted to this condition. 

3.2. Numerical methods 

Gas Evolution: The general relativistic Euler system of 
equations (5), is solved in time using the method of lines, 
that uses a third order total variation diminishing (TVD) 
Runge-Kutta time integrator (Shu & Osher 1988). These 
are discretized using a finite volume approximation to¬ 
gether with high resolution shock capturing methods. 
Specifically, we use the HLLE (Harten et al. 1983; Ein- 
feldt 1988) approximate Riemann solvers in combination 
with the minmod linear piecewise reconstructor. The nu¬ 
merical fluxes and sources in (5) depend both on the con¬ 
servative and on the primitive variables w = {po^v\p). 
Then, in order to express primitive in terms of conserva¬ 
tive variables, we use a Newton-Raphson algorithm each 
time step within the evolution scheme. 

Boundary Conditions: We study numerically the 
relativistic gas on the equatorial plane, in the do¬ 
main [rexc^ ^max] X [0,27r) with resolution (Ar, A0) = 
(0.158,0.05) for all cases. We choose the interior bound¬ 
ary rexc to be inside the black hole horizon, were we 
apply a numerical excision (Seidel & Suen 1992), i.e., we 
apply a cutoff inside the event horizon which is possible 
due to the Kerr-Schild coordinates used by CAFE. The 
exterior boundary rmax^ is splitted in to two halves, one 
in which the gas enters the domain where we apply in¬ 
flow boundary conditions, and a second half where the 
gas leaves the domain and we apply outflow boundary 
conditions there. Besides, in the angular domain we use 
periodic boundary conditions. 

Initial data: As initial data, we consider a wind, that 
fills the whole domain, moving on the equatorial plane 
along the x direction with non constant density profile. 
We characterize the initial velocity field v'^ in terms of 
the asymptotic initial velocity Voo as done in Cruz-Osorio 
et al. ( 2012 ); Font et al. (1999), where the relation v‘^ = 


Viv'^ = is satisfied. Using Kerr-Schild coordinates the 
explicit expressions for the velocity vector field v'^ are 
given by: 


= HiVooCoscI) + H2VooSm(l), ( 10 ) 

= —H^VooSincI) H^VooCoscf). ( 11 ) 

where Hi {i = 1 , 2 ,3,4) represent functions associated 
with the components of the three metric. 


H^ = 
H2 = 
Hs 

H. 


's/^rr 


H\'^rr T H^'^j^(j) 


H\^rr T -^^477^0 


("^1T -^^4 700 T 277iiE47r0) 


( 12 ) 

(13) 

^4) 


27r0 


(15) 


The profiles of the components of the field velocity can 
be seen in figure 1 . 

Following Ruffert’s work (Ruffert & Anzer 1995; Ruf- 
fert 1997, 1999), the initial density gradient is chosen in 
such way that it is perpendicular to the gas motion and 
also proposed as an hyperbolic function in order to serve 
as a cutoff at large distances for large density gradients. 
The density distribution, in polar coordinates, is given 
by the following expresion 


Pini = Pojl - ^tanh 


, (rsin [7 +acos[ 7 ) i-| 

a 


where po is a constant density value, a = J/M is the 
spin of the black hole, Cp is the parameter specifying the 
magnitude of the density gradient and is the accretion 
radius, which is defined in terms of the asymptotic values 
of the sound speed Cg^ as Petrich et al. (1989) 


Ta = 



(17) 


Figure 2 illustrates the density profile for two values of 

Once we fix the value of Cgoo and assume the gradient 
density profile (16), the pressure can be found from the 
equation of state as: pini = c^poPimlil' - where 

7 i = 7 /( 7 ~l)- III order to avoid negative and zero values 
on the pressure, the condition Cgoo < ^7 ~ I fo be 
satisfied. 

Now, Voo and Cgoo are two useful parameters which 
define the relativistic Mach number at infinity, = 

r^’oc/TsCsoo = TAdoo/Ts. Here Fg is the Lorentz 
factor calculated with the speed of sound and A4oo 
is the asymptotic Newtonian Mach number used to 
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Fig. 2. — This figure shows the initial density profile for two values of the density gradient parameter, Cp = 0.2 and Cp = 0.5. As can be 
observed, for higher values of Cp the drop in density is considerably more pronounced. We remind the reader that length is in units of M 
with 1 Mq = 1.48 X 10^ cm. 

















7 


parametrize the initial configurations. The initial data is 
parametrized using the Mach number given in Table 1. 

We use CAFE (Lora-Clavijo et ah 2015), a fully three- 
dimensional relativistic-magnetohydrodynamic (MHD) 
code. Although the MHD in CAFE is written for a 
Minkowsky space-time (it does not use curved space- 
time), the hydro dynamical (HD) routine does utilize 
the equatorial (Cruz-Osorio et al. 2012, 2013) and axial 
(Lora-Clavijo & Guzman 2013) symmetries in order to 
allow simulations using fixed curved space-times. CAFE 
also solves the Einstein field equations coupled with rel¬ 
ativistic HD in spherical symmetry (Guzman et al. 2012; 
Lora-Glavijo et al. 2013b, 2014b). The numerical meth¬ 
ods available in GAFF include several reconstructors. 
For instance MINMOD and MG linear piecewise meth¬ 
ods; for higher reconstructors GAFE uses PPM parabolic 
method and WEN05 polynomial method. All of these 
reconstructors are used in combination with HLLE ap¬ 
proximate Riemann solver flux formula. Goncerning rel¬ 
ativistic MHD, and in order to preserve magnetic field di¬ 
vergence, GAFE uses flux constraint-transport and diver¬ 
gence cleaning methods. Numerical tests of the numeri¬ 
cal implementation in GAFE can be found, for the rela¬ 
tivistic HD, in Lora-Glavijo et al. (2013a); Lora-Glavijo & 
Guzman (2013) and for relativistic MHD in Lora-Glavijo 
et al. (2015). 


Aloo 

a=0 

ep Ta 

a=0.3 

ep Ta 

a=0.5 

ep Ta 

a=0.7 

ep Ta 

a=0.9 

ep Ta 

5 

0.0 

0.2 3.85 

0.5 

0.0 

0.2 3.85 

0.5 

0.0 

0.2 3.85 

0.5 

0.0 

0.2 3.85 

0.5 

0.0 

0.2 3.85 

0.5 

4 

0.0 

0.2 5.88 

0.5 

0.0 

0.2 5.88 

0.5 

0.0 

0.2 5.88 

0.5 

0.0 

0.2 5.88 

0.5 

0.0 

0.2 5.88 

0.5 

3 

0.0 

0.2 10.0 
0.5 

0.0 

0.2 10.0 
0.5 

0.0 

0.2 10.0 
0.5 

0.0 

0.2 10.0 
0.5 

0.0 

0.2 10.0 
0.5 

2 

0.0 

0.2 25.0 

0.5 

0.0 

0.2 25.0 

0.5 

0.0 

0.2 25.0 

0.5 

0.0 

0.2 25.0 

0.5 

0.0 

0.2 25.0 

0.5 


TABLE 1 

In this table, we summarize the models studied in figures 3 
AND 4. All the different configurations considered here 

ASUMME THE RELATIVISTIC SOUND SPEED AT INFINTY IS Cgoo =0.1 
AND ADIABATIC INDEX 7 = 5/3. 


3.3. Diagnostics 

In order to diagnose the amount of mass and angular 
momentum accreted by the Kerr BH we implement a 
detector located as close as possible to the (outer) event 
horizon at: 

r+=M+ yM2 - a2. 

This implies we define a sphere where we compute said 
scalars. The following formulae represent both quantities 
respectively: 

p27T 

M= / ay/^D{v'^ — 13'^/a)d(t), (18) 

Jo 

p27T nVdet 

P^ = - ay/^T^^d(t)+ / / S'*‘drd(t), (19) 

Jo Jo Jrexc 


where stands for the (j) component {j = (j)) of the 
source term of equation (8), Vexc is the excision radius 
and Vdet is the radius where the detector lays. 

4. RESULTS 

We have studied the parameter space for mass and 
angular-momentum accretion rates. Our parameters are 
the Mach number. Ad (the columns in our plot array in 
figures 3 and 4), the spin of the BH, a (the rows in our 
plot arrays in figures 3 and 4) and the density-gradient 
parameter of the gas, Cp (represented by the different 
color lines in our plot array in figures 3 and 4). 

All the shown figures (5) are snapshots of the system 
in steady state, or the stationary regime. These figures 
illustrate the BH rotating counter-clockwise for positive 
values of a and with the flow moving from left to right. 
That is, they correspond to timescales where mass and 
angular momentum accretion rates have stabilized and 
thus the curves (in figures 3 and 4) have plateaued. 

Our results, as seen in the first column (for Cp = 0) of 
figures 5 correctly reproduce the expected properties of 
the shock cone when there exists no density gradient in 
the medium where the black hole moves. I.e., we observe 
the formation of a symmetric shock cone whose width 
depends on the Mach number; as the Mach number in¬ 
creases the shock-cone angle decreases (Font & Ibanez 
1998a; Font et al. 1998; Font & Ibanez 1998b; Font et al. 
1999; Gruz-Osorio et al. 2012; Lora-Clavijo & Guzman 
2013). 

4.1. Morphology 

We present, for the first time, the 2D morphology of 
the relativistic BHL accretion onto a BH considering den¬ 
sity gradients (figures 5). It is worth mentioning that in 
order to illustrate the general morphology of the system, 
we present only the case of Al = 5; however we have cov¬ 
ered all configurations prented in this paper. This first 
attempt is done in slab symmetry, which is the first of 
a series of steps towards more realistic 3D simmulations. 
The color gradient represents the logarithm of the gas 
density in geometrical units normalized to the BH mass. 
The contour lines simply emphasize the density gradient. 
These figures dramatically illustrate the effect that the 
density gradient has on the shock cone once steady state 
is acheived. 

From the first column of figures 5 we observe that (as 
in Font et al. 1999, where zero density gradient is consid¬ 
ered) as a increases so does the induced angular momen¬ 
tum in the wind. Now, as the density gradient increases, 
we observe the most notable feature in these figures, i.e., 
the shock cone is further pushed towards the lower den¬ 
sity region. 

4.2. Mass and Angular Momentum Accretion Rates 

From the plots on figures 3 and 4, we observe that, as 
expected, the mass accretion rate decreases as the Mach 
number, Aloo, increases. Unlike the Newtonian case 
where no trend between the mass-accretion-rate fluctu¬ 
ations and the density gradient was discerned (Ruffert 
1999) our results exhibit a stationary state rather quickly 
and show that as Cp increases the mass accretion rate 
slightly decrases. The faster the fluid moves with re¬ 
spect to the BH the faster the steady states settle in. 















The mass accretion rate is, mostly, independent of the 
spin parameter of the BH. The reason for the the steady 
state settling is rather quickly may be due to the fact 
that our simmulations are done in slab symmetry as op¬ 
posed to 3D. It could also be that we are dealing with 
relativistic flows. We confirm that the accretion rates 
do decrease slightly when the density gradient increases. 
Besides, when the BH spin increases the effects of the 
wind density gradient on the mass accretion rate become 
slightly stronger. 

We further explore cases where the accretor has an¬ 
gular momentum. As Aloo decreases the difference be¬ 
tween the mass accretion rate, for different e^, becomes 
larger. This is not noticeable, however, in the case where 
AIoo = 2, because for lower velocities steady state settles 
in much later (notice the timescale is ^ 3 times longer 
and it is still not fully achieved; this would be consistent 
with the Newtonian results in Ruffert 1999). 

The behavior of the angular momentum accretion rate 
(figs. 4), instead, clearly shows that for higher wind ve¬ 
locities (increasing Moo) steady state is reached much 
more quickly. For Moo = ^ and Moo = ^ steady state is 
achieved after about t = 200 and t = 100, respectively, 
regardless of BH spin or density gradient. Instead, look¬ 
ing at Aloo = 3 and, espacially at Atoo = 2 it is clear that 
a steady state is not fully reached in t = 600 for the for¬ 
mer and not even in t = 2000 for the later. This trend is 
more noticeable for lower density gradient and/or higher 
BH spin. In other words, it appears that increasing the 
density gradient helps stabilize the angular momentum 
accretion rate at an earlier time. We also observe that 
the angular momentum rates extered on the wind de¬ 
crease into negative values further and further the more 
the BH spin, a, increases. It can also be seen that higher 
Mach number has the effect of slightly increasing the an¬ 
gular momentum accretion rate at steady state for large 
density gradient (e^ = 0.5). For lower density gradient 
(ep = 0.2) this trend is not evident. And, for no density 
gradient (e^ = 0) there seems to exist no correlation. 

We have produced, as well, a set of simulations where 
the spin of the BH counterrotates, a < 0. The first case 
has a = —0.9, 7 = 5/3 and AIoo = 0-3; the second case 
has a = —0.9, 7 = 5/3 and Aloo = 0-4; and, finally the 
third case has a = —0.99, 7 = 5/3 and Aloo = 0-4. In all 
three cases we performed simulations for the same three 
values of the density gradient, i.e., Cp = 0, 0.2 and 0.5. 
Opposite to the cases described above, here the angular 
momentum of the accreted material is parallel to that of 
the BH, thus, it is expected that mass accretion will tend 
to increase the spin of the BH on the long run. 

As can be observed from the simulations for low Mach 
number (see fig. 7), a disk-like structure forms. The ve¬ 
locity field has positive radial coordinate, thus, little or 
no accretion occurs. 

5. DISCUSSION AND CONCLUSIONS 

We have performed a parameter study of BHL accre¬ 
tion of a relativistic wind with density gradient onto 
Schwarschild and Kerr BHs with different spin param¬ 
eters. A discussion of our results follows in the next 
paragraphs. 

Comparing figures 3 vs 4, and 6, it is interesting to 
note that the mass accretion rate, unlike the angular mo¬ 
mentum rate (see first column plots in fig. 4), does not 


depend strongly on the wind density gradient regard¬ 
less of Mach number or BH spin. As can be observed 
from our fig. 3, observing the plots from left to right, 
as the Mach number increases the mass accretion rate 
decreases; in agreement with previous studies as well as 
theory. It is also important to note that, as the Mach 
number increases, steady state is obtained much more 
rapidly. This is further implied by the results in Ruffert 
(1999), where only Newtonian velocities are acheived and 
the plots show much more variabilty in both mass and 
angular momentum accretion. 

From the plots in figure 4 we find that the wind may be 
accreted with low or even negative values of angular mo¬ 
mentum with respect to the spin of the BH. This implies 
that, were a significant amount of material to be accreted 
(over long periods of time), the spin of the BH could be 
brought down significantly or even reversed in the cases 
where the angular momentum has negative values. On 
the other end of the BH spin spectrum, from fig. 6, it is 
clear that the spin of a BH can increase over time if the 
spin and the wind density gradient are properly aligned. 

From the morphology plots (figs. 5), we observe a new 
signature, i.e., the Mach cones wrap around the BH, even 
with low or nil spin, due to the density gradient of the 
winds. 

Probably, due to the fact that our runs do stabilize 
rather quickly, we observe that an accretion disk starts 
to form for 7 = 5/3, low Mach number, high BH spin 
and high density gradient (e^ = 0.5; see Fig. 7). 

As noted above, in the simulations where a disk forms 
the velocity field has positive radial component, which 
would imply no accretion occurs. However, the fluid 
modeled by our code has no account of radiative losses. 
Thus, in a more realistic simulation it is highly likely that 
energy would be radiated away and angular momentum 
would be transported outwards within the disk allowing 
for substantial accretion onto the BH. 

5.1. Astrophysical Scenarios and Applications 

From figure 3 we can observe that for a ^ IOMq 
BH steady-state accretion can be acheived after a few 
milliseconds. For IMBHs the timescales necessary for 
reaching steady-state accretion go from seconds to hours. 
Whereas for SMBHs said timescales go from hours to 
years. 

In the hypernova-explosion scenario for sBHs and 
IMBHs, accretion will probably last a few minutes thus 
little mass or angular momentum may be accreted. 

For the scenario where a SMBH crosses through the 
jet of a companion SMBH substantial mass may be ac¬ 
creted if the orbit and the supply of material to the BH 
producing the jet are stable. However, the effects of the 
density gradient modifying the spin of the BH may be 
cancelled out if the axis of the jet is on the orbital plane. 
This occurs because as the BH enters the cone of the jet 
the density increases, however, as the BH exits the jet, 
the opposite effect takes place, thus cancelling most of 
the effect. If the axis of the jet is slightly off the orbital 
plane then there is a component of density gradient that 
does not get cancelled out and a BH spin can be built 
up over many orbital periods. If two jets, as opposed 
to one, exist and they are symmetric, the effect of the 
density gradient on the BH spin will also be cancelled 
out. Were the jet to process (so that the axis of the jet is 
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sometimes above the orbital plane and sometimes below) 
on a timescale much larger than the orbital period spin 
reversal could be observed on the BH that cuts through 
the jet. Maybe more important is the fact that these 
processes are accreting mass but more likely than not, 
no angular momentum, thus the spin, a, of the BH will 
decrease (as the spin a oc J/M‘^). 

The reader may ask how likely it is that the jet of a 
SMBH may hit another one. There are a few examples 
where binary SMBHs have been observed (e.g., Valtonen 
et al. 2008; Fabbiano et al. 2011; Valtonen et al. 2012; 
Graham et al. 2015). The most relevant parameters to 
estimate the likelihood of such an event, we suspect, are 
the following: First, the jet must have a large angle, oth¬ 
erwise chances are really small. Second, it is suspected 
that these binaries are the product of a merger of two 
galaxies, hence it is likely that if the spin of the SMBHs 
are similar to those of their hosts (e.g., Barausse 2012) 
and they collide at random angles, the SMBHs spins will 
also be random. Third, the impact parameter of the colli¬ 
sion should, preferably, be small such that the BHs do not 
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have to travel far within the newly merged galaxy, oth¬ 
erwise they may accrete substantial amounts of matter 
(see, e.g., Armitage & Natarajan 2002; Dotti et al. 2009) 
which will have a prefered angular momentum (that of 
the host galaxy) and may reorient the axes towards that 
of the reshaped host (Dotti et al. 2009), and, thus, pre¬ 
venting them from hitting each other with their jets. The 
mass of the SMBHs will be important as well, as the spin 
of a more massive BH will be less affected by accretion 
as it drifts inwards to meet the companion BH (Barausse 
2012). 
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Fig. 3.— This figure shows the evolution of the mass accretion rates (the mass accretion rate is rescaled by 10“^®). The mass accretion 
rate is measured near the event horizon of the BH (whose radius is a function of the spin parameter). Each model is run, at least, until 
steady state is acheived. We show the cases considering various values of the Mach number A4 = 2, 3, 4, 5 (columns) and different values 
for the BH spin parameter a = 0.3, 0.5, 0.7, 0.9 (rows). In each figure we consider three values of the density gradient Cp = 0, 0.2, 0.5. 
We find that the mass accretion rate increases as the Mach number decreases. In contrast, we can see that there is litle infiuence of the 
density gradient on the accretion rate. We remind the reader that time is in units of M with I Mq = 4.93 x I0“® s. 
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Fig. 4. — In this figure we show the angular momentum accretion rates vs time (the angular momentun rate is rescaled by 10“^^). Again 
we we show the models corresponding to Mach numbers A1 = 2, 3, 4, 5 (columns) and for angular parameter a = 0.3, 0.5, 0.7, 0.9 (rows). 
For each value of the Mach number we study three values of the density gradient Cp = 0, 0.2, 0.5. We can observe notable changes in the 
angular momentum rates when the density gradient parameter increase. All rates are measured near the event horizon. We remind the 
reader that time is in units of M with 1 Mq = 4.93 X 10“® s. 
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Cn = 0 


= 0.2 



Fig. 5. — Morphology of the rest mass density at stationary state, we can observe the shock cones are dragged due to the density 
gradients. We show models with Mach number A4 = 5. Different rows, from top to bottom show results for a = 0, 0.5, 0.9. The columns, 
from left to right show density gradient parameter Cp = 0, 0.2, 0.5. The contour lines range is [—12.6, —10] and the space between a line 
and another is 0.1 in each plot. We remind the reader that length is in units of M with 1 Mq = 1.48 x 10^ cm. 
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AIqo = 3, a = —0.9 


AIqo = 4:, a = —0.9 


M-oo = 4:, a = —0.99 





Fig. 6. — In this figure we show the angular momentum rates of the special cases considering counterrotating black holes and density 
gradients Cp = 0, 0.2, 0.5. In contrast with positive values of the angular momentum of the black hole we can observe that, in all cases, 
the angular momentum rates are positive. 



l 





t t 

Fig. 7. — We show, on the top two figures, the morphology of the rest mass density and the velocity field vectors. We show two close-ups, 
one at r = 8 (left) and other at r = 20 (right). On bottom we show the evolution of mass (left) and angular momentum (right) accretion 
rates. This special case corresponds to a model with paramenters 7 = 5/3, A4oo = I, a = 0.5 and Cp = 0.5. In this model we find the 
formation of a disk-like structure (see the 2D figures in the upper panels), the velocity-field vectors shows cuasi-circular trajectories of of 
the fiuid elements. Notice that steady state has not been fully acheived for these cases after 4,500 M. 
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